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ABSTRACT 

Pulsar timing arrays (PTAs) are expected to detect gravitational waves (GWs) from 
individual low-redshift (z ^ 1.5) compact supermassive {M > lO^M©) black hole 
(SMBH) binaries with orbital periods of ~ 0.1 — 10 yr. Identifying the electromagnetic 
(EM) counterparts of these sources would provide confirmation of putative direct de- 
tections of GWs, present a rare opportunity to study the environments of compact 
SMBH binaries, and could enable the use of these sources as standard sirens for cos- 
mology. Here we consider the feasibility of such an EM identification. We show that 
because the host galaxies of resolved PTA sources are expected to be exceptionally 
massive and rare, it should be possible to find unique hosts of resolved sources out 
to z « 0.2. At higher redshifts, the PTA error boxes are larger, and may contain as 
many as ~ 100 massive-galaxy interlopers. The number of candidates, however, re- 
mains tractable for follow-up searches in upcoming wide-field EM surveys. We develop 
a toy model to characterize the dynamics and the thermal emission from a geometri- 
cally thin gaseous disc accreting onto a PTA-source SMBH binary. Our model predicts 
that at optical and infrared frequencies, the source should appear similar to a typi- 
cal luminous active galactic nucleus (AGN). However, owing to the evacuation of the 
accretion flow by the binary's tidal torques, the source might have an unusually low 
soft X-ray luminosity and weak UV and broad optical emission lines, as compared to 
an AGN powered by a single SMBH with the same total mass. For sources at z ^ 1, 
the decrement in the rest-frame UV should be observable as an extremely red optical 
color. These properties would make the PTA sources stand out among optically lumi- 
nous AGN, and could allow their unique identification. Our results also suggest that 
accreting compact SMBH binaries may be included among the observed population 
of optically bright, X-ray-dim AGN. 

Key v^rords: black hole physics — gravitational waves — accretion, accretion discs 
— galaxies: active 



1 INTRODUCTION 

Over tlie last several years, the possibility of ob- 
serving both the gravitational-wave (GW) and electro- 
magnetic (EM) emission signatures of coalescing super 
massive black hole (SMBH) binaries has received 



tense attention 
I2OO7I : iKocsis. Haiman 



JHolz fc Hughed I2OO5I; iKocsis et all 



Menou 



iDotti et all 



200e, 



200e 



for specific proposed mecha nisms for EM signatures, see 
Armitage fc NataraianI |2002| and iMil osavlic vic fc PhinnevI 

• iHai n 



2OO5I. as well as recent reviews by iHai man c t al.l |2009| and 
Schnittmanll201ll ). The bursts of GWs emitted by such sys- 



tems can now be predicted by numerical general relativity 



E-niail:taka@mpa-garching. mpg.de 



(|Pretoriusll2"005l : iBaker et al1l2006l : ICampanelli et al.ll2006l ). 
and are expected to be observed by current and future detec- 
tors. The temporal evolution of the gravitational waveform 
can be used to extract the luminosity distance, help con- 
strain the location of the source on the sky, and determine 
the masses and spins of the SMBHs. If an EM signature 
of the coalescence can also be identified, this would allow 
for a determination of the source redshift, turning merg- 
ing black holes into "standard sirens" for probing cosmic 
expansion|j Such multi-messenger observations would also 
enable astronomical investigations of SMBHs whose masses, 



The importance of su ch GW+EM o bservations for cosmogra- 
phy was first discussed bv lSchutj 1119861) in the context of merging 
neutron star binaries. 
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spins and orbital parameters are already known, presenting 
ideal laboratories for investigating accretion physics in ac- 
tive galactic nuclei (AGN). Furthermore, if m ajor mergers of 



galax i es trigger luminous AGN activity (e.g., Sanders et al 



1988| : lHernauistlll989l: ICarlberd [iQQOI: iBarnes fc Hernguisl 



199ll: iHernguist fc MihosI Il995l: iMihos fc Hernguistl llQol 



Kauffmann fc Haehneltl I2OO0I : iHopkins et al.l |2007| . i20o'i ). 
then the characteristic EM emission promptly following 
the SMBH coalescence niay her ald the birth of a quasar 
iTanaka. Haiman fc Menouflioiol ') . 

To date, theoretical studies of EM signatures of GW- 
emitting SMBH binaries have largely centred on systems 
predicted to be detectable by future space-based laser in- 
terferometers such as Laser Interferometer Space Antennae 
{LISA), with total mass ~ 10^"'^(1 + z)~^Mq and redshifts 
as high as 2 > 10. For the purposes of multi-messenger 
astronomy, a paramount feature of interferometer-detectors 
is the precision with they c ould determine the angular sky 
position of SMBH sources (IKocsis et al.l 120061 1: to ^ 1 deg 
(|Vecchid |2004| : iLang fc Hughes! |200~ ), or perha ps even to 
^ 1' when spin- induced pr ecession (Lang fc Hughe s 200|j) 
or higher-order harmonics (|Mc Williams et al.ll2010l ) are in- 
cluded in the analysis of the waveform. 

In this paper, we evaluate the prospects of perform- 
ing multi-messenger using pulsar timing arrays (PTAs), 
which aim to detect low-frequency (nHz) GWs through 
precision timing of Galactic millisecond pulsars. A ma- 
jor goal for PTAs is to detect the stochastic GW back- 
ground due to the population of comp act SMBH b i naries 
in our cosmic neighborhood {z ^ 2). Ijenet et al.l (|2005l ) 
showed that this requires the timing of some > 20 pulsars 
to ~ 100 ns precision over 5 years, a g oal that could be 
achieved by currentl y operational arrays (jManchesteiiuOOg : 
IVerbiest et al.ll2009l . and refs. therein). Recent theoretical 
popu l ation-synthesis studies (ISesana. Vecchio fc Volonteril 
I2OO9I : ISesanafc Vecchidl2010l : iKocsis fc Sesanall201ll ') haw 
predicted that PTA observations may also be able to in- 
dividually resolve the most massive and/or most nearby 
SMBH binaries that stick out above the stochastic back- 
ground|j The primary factor that determines the number 
of detectable individual sources is expected to be the array 
sensitivity, with theoretical uncertainties such as the binary 
mass function and orbital evolution affecting predictio ns of 
the o bservable population by factors of unity ( Scsana c t al.l 

The systems individually detected by PTAs will differ 
from LISA targets in several ways. First, according to the 
population synthesis studies, the binaries will have chirp 
masses 

M = mI^'^mI''"M-^'^ = rf^^'M, (1) 

typically around M ~ 10*'^ M©. Above, M2 < Mi are the 
masses of each member of the binary, and 77 = MiM2/M'^ < 

■^ In principle, PTAs can also detect GW-nie mory bursts 
from individual major SMBH mergers jSetol |2009|: 



Pshirkov. Baskaran fc Postnovl I2OICII ; Ivan Haasteren fc LevinI 
201Cri . However, the number of detections is expected to be 



quite low, perhaps much lower th an unity ov er a full observation 
campaign (see last paragraph in ISetd 12009 ) . We therefore limit 
the present paper to the population and EM signatures of 
compact SMBH binaries. 



1/4 is the symmetric mass ratio. In other words, the systems 
of interest here are much more massive than those relevant 
for LISA. Second, PTA sources lie at much closer cosmolog- 
ical distances. While the redshift probability distribution is 
poorly known, it is expected to decline steeply outside the 
range 0.1 ^ z ^ 1.5. The low cutoff is due to the smaller 
volume contained within z, whereas the high cutoff is due 
to the attenuation of the GW signal and the decline of the 
intrinsic SMBH merger rate. Third, the PTA targets are not 
as compact, having periods of P ~ 1 yr. The larger separa- 
tions correspond to slower orbital decay; the vast majority 
of PTA sources will not coalesce within a human lifetime. 
That the orbit — and thus the waveform — is expected 
to evolve much less appreciably with time poses a particu- 
larly difhcult challenge for determining the masses and lu- 
minosity distances of PTA sources. Particularly important 
for the prospects of performing synergistic EM observations 
on these sources is their relatively poor sky-localization. The 
solid-angle uncertainty in the source position may be any- 
where from AQ, ^ 3 deg^ (in the best-case scenario in which 
the co ntributions to the signal from individual pulsars are 
known: ICorbin fc Cornishll20ld . hereafter CCIO) to as large 
as AQ, ~ 40 deg'^ (if individu al pulsar contributions c annot 
be extracted from the data; ISesana fc Vecchid l20ld . here- 
after SVIO). 

Several tell-tale EM features of such compact SMBH 
binaries have been proposed in the literature, including 
periodic emission modula ted at the orbital frequency 
of the PTA source (e.g., iHaiman et al.l |2009| and refer- 
ences therein) and double-peak e d br o ad emission lines 
(e.g. iGaskelJ Il99d: IZhou et al.l |2004|: iBogdanovic et all 



200i 



BlechafcLoed I2OO8I: iBoroson fc Laueij |2009|: 



Bogdanovic. Eracleous fc Sigurdsson 20091: iDotti et al.l 



20091 : iTang fc Grindlavl I2OO9I : IShen fc Loed l20ld ). In 
this work, we will exp lo re in detail the suggestion by 
iMilosavlievic fc PhinnevI (|200a ) that tidal torques of a 
near-merger SMBH binary will have evacuated the central 
region of its accretion disc, resulting in a markedly softer 
thermal spectrum. 

Below, we investigate whether individually resolved 
PTA sources may be viable targets for EM identification. 
In particular, we address the following two questions: 

(i) What is the average number Ng of candidate host 
galaxies — that is, interloping galaxies that could plausibly 
harbour the GW source — in a typical error box of a PTA 
detection? Of particular interest is whether there are plausi- 
ble scenarios for detecting individually resolved sources with 
Ng < 1, i.e., cases where the source may be uniquely identi- 
fied with an EM search of the three-dimensional PTA error 
box. 

(ii) In cases where A*'g > 1, what can be done to dis- 
tinguish the true host galaxy of the source from the other 
interlopers? Motivated by the hypothesis that galaxy merg- 
ers can fuel AGN activity, we will consider the differences 
in thermal emission properties predicted by disc models of 
AGN powered by a compact SMBH binary as opposed to 
one powered by a solitary SMBH of the same total mass. 

This paper is organized as follows. In ^ we provide a 
brief overview of the expected population of PTA-resolved 
SMBH binaries as well as the anticipated detection error box 
of such objects. We consider the error box prescriptions of 
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ISVlCl and ICClOl . characterize the types of astronomical ob- 
jects that are plausible hosts of a PTA-resoIved binary, and 
estimate the number of such objects. In §3, we describe a toy 
model to calculate the dynamical state and thermal emis- 
sion features of gas accreting onto a resolved SMBH binary. 
We discuss several features predicted by the model which, 
if observed, could help the EM identification of individually 
resolved PTA sources. We summarize our findings and offer 
our conclusions in §4. 

Throughout this paper, c denotes the speed of light; G 
is the gravitational constant; fee is the Boltzmann constant; 
asB is the Stefan-Boltzmann constant; and m^ is the mass 
of the proton. 



2 PLAUSIBLE HOSTS OF PTA-RESOLVED 
BINARIES 

As stated in iJT] theoretical models predict that SMBH bi- 
naries individually resolved by PTAs are most likely to have 
masses of M > 1O*M0, observed periods of Pobs ~ 1 yr, and 
lie in a redshift range 0.1 ;$ z ;$ 1.5. Given the PTA detec- 
tion of such a GW-source binary, we wish to evaluate the 
number Ng of candidate galaxies that could plausibly host 
it. To this end, we will first review the volume of the error 
box in which we must look for the source, based on previous 
work on the source localization capability of PTAs. Then, 
we will evaluate the number of interloping host galaxies in 
the error box by estimating the number of (i) sufficiently 
massive dark matter halos, (ii) sufficiently luminous lumi- 
nous galaxies, and (iii) AGN. We adopt a standard ACDM 
cosmology with dimensionless Hubble parameter h = 0.70, 
matter density Jim = 0.27, dark energy density Q.a ~ 1 — fim, 
baryon density ilh = 0.046, and p ower spectrum amp litude 
as = 0.81. (WMAP 7-year resuHs. [jaxosik et al.|[201ll '). 



2.1 The PTA Error Box 

We consider two different estimates of the size of the er- 
ror box of PT A-resolv ed sources. The first is based on the 
calculations bv lSVld . who assumed that the contributions 
to the signal from GW perturbations at the individual pul- 
sars (the so-called "pulsar term") cannot be extracted from 
the PTA data. Those authors found that for a hypothetical 
array containing 100 pulsars and covering the whole sky, a 
typical resolved source could be be localized within a solid- 
angle error of A^ ~ 40 (SNR/10)"^ deg^, where SNR is the 
signal-to-noise ratio. The fractional error of the signal am- 
plitude A ex M^'^Dl^ wiU be of order ~ 30 (SNR/10)"^%. 
These values are statistical means from their Monte-Carlo 
simulations; the localization and amplitude measurement for 
an individual source may be better or worse, depending on 
the specific parameters of the binary and the orientation of 
the pulsars with respect to the source. Given the wide spread 
in chirp mass distribution predicted by population synthesis 
models of resolved sources, in the absence of an independent 
measurement of M the only constraint on Dl comes from 
the maximum distance at which PTAs are expected to de- 
tect ind ividually resolved sources. The population synthesis 
mode ls (|Sesana. Vecchio fc Volonterl2009l : lKocsis fc Sesanal 
I2OIII ) predict that the majority of resolved sources will lie 



below a maximum redshift Zmax ~ 1-5, or a luminosity dis- 
tance below Dl, max ^ 10* Mpc. This "worst-case" error box 
has a comoving volume of 



Al/ 



CSV) 



3 X lO'' 



Afi 



Mpc'^ 



(2) 



,40deg2 

More optimistic numbers are obtained by ICCld . who 
suggested that utilizing information on the distances to in- 
dividual pulsars in the array can greatly enhance the mea- 
surement capabilities of PTAs. They concluded that if the 
individual pulsar term can be extracted from the signal, 
then this would double the signal power and enable di- 
rect measurement of the chirp mass. They estimate that 
for a system with SNR= 20 (corresponding to a detec- 
tion of SNR= 10 without pulsar distance information), a 
resolved source can be localized with distance and angu- 
lar errors of /S.Dl/Dl < 20% and AfJ < 3deg^, respec- 
tively. Noting that the comoving distance D{z) in the rel- 
evant redshift range can be analytically approximateqj as 
D ~ c Hn z (1 — 0.22), we may estimate the error box in the 
ICCld scenario as 



AF'*^*^^ « 1.2x10^ 2^(1-0.22)-^ 



2.2 Interloper Counts 



3degV V 20% J ' 

(3) 



Where does an individually resolved PTA source live? 
The mass M of a nuclear SMBH is known to cor- 
relate with the velocity dis persion a of the host 
galaxy (the "M - a r elation" : iFerrarese fc Merritd l200d : 
iGebhardt et al.l l200d : iTremaine et al.1 I2OO2I ). as well as 
with t he stellar luminosity of the host (the "M — L rela- 
tion" : iKorm endv fc Richstondll995l:lMagorrian et al.lll998l : 
iHaring fc R ix 200^^; iLauer et al.l 120071 ): with more massive 
halos and luminous galaxies hosting more massive SMBHs. 
That resolved PTA sources are expected to be exceptionally 
massive {M > Kf Mq) implies that the host should be a 
giant elliptical galaxy or be among the most massive spi- 
ral galaxies (with veloci ty dispersion a > 200 km s"'^ of the 
spheroid component; e.g lGiiltekin et al.ll200d ). 

Assuming that SMBH b inaries are able to overcome the 
"fina l parsec" problem (e.g..lEscala et al.ll2005l:lMaver et al.l 
l2007l : ICallegari et al.|[200^ 'Colpi et al.ll200d : lHavasakill2009r 
see, however. iLodato et al ., 20091. we expect the PTA host 
galaxy to be the product of a relatively recent merger. A 
natural question to ask is whether such galaxies typically 
lie in the field, or in the centre of a cluster. We can answer 
this question qualitatively by considering the dependence of 
the major merger rates of the most mass ive dark matter 
halos on t heir environments. Analyses by iFakhouri fc Mai 
(|2009l ) andlBonoh et al.l (l20ld) of the MiUennium simulation 
results (|Springel et al.ll2005l ) indicate that while the rate of 
major mergers is enhanced in over-dense environments, this 
effect is weak: for halo masses and redshifts of interest [M > 
10^^ Mq and at 2 SS 1-5), the ratio of merger rates between 

^ This fitting formula has an error of less than 1% in D at 
2 < 1.4 and roughly 5% at z = 1.9. It is provided for the reader's 
convenience; all distance and volume calculations in this paper 
are performed using exact expressions. 
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the most and least over-dense regions is of order unity. We 
interpret this result to mean that there is no strong reason 
to search for PTA sources in galaxy clusters as opposed to 
those in the field. 



2.2.1 The most massive halos 

One conservative way to estimate the number Ng of host 
galaxy candidates in the error box is to simply count the 
dark matter halos that are massive enough to plausibly har- 
bour th e source SMBH bina ry. We use the observational re- 
sults of iDutton et al.l (|20ld ). who infer a double-power- law 
fit for the relation between the SMBH mass M and the host 
halo mass Afhaio for local elliptical galaxies. We extrapolate 
their results to higher redshifts by postulating the canonical 
z-dependence b ased on the theory for spherically collapsing 
halos (see, e.g.. lWvithe fc Loebl[20oi ). 



Mhalo(M,z)cxF(2) 



d{z) Ae(0) 
d(0) Ae(^)' 



(4) 



where d(2) = -[(Qm/nA)(l + z)^ + l]"^ and Ac(z) = ISyr^-f 
82d(z) - 39d^{z). We obtain 



Mhaio « 2.3 X lO'^'M^'-^'' 



13Mq 



F{z) Mq, (5) 



where Mg = M/(10®Mq). We estimate the number of can- 
didate host halos inside the three-dimensio nal PTA error 
box b y integrating the halo mass function of I Jenkins et al.l 
(|200ll : their equation 9) above A^fhaio- 

The most massive halos with Mhaio ^ few x 10^'* Mq, 
which are associated with galaxy clusters, may be ex- 
pected to contain more than one plausible host galaxy. 
Since the halo mass function at Mhaio ^ few x 10^'* Af© 
drops much more steeply than linear with mass, whereas 
the sub-halo mass func tion increases less steeply than lin- 
ear (jGiocoli et al.l [2010|) , most galaxies with halos masses 
~ lO^^Af© will reside in the field, rather than in groups and 
clusters. The multiple occupancy of massive galaxies in the 
most massive halos will then represent only a small increase 
in our total counts of interlopers. As the purpose of the ex- 
ercise in this section is to give order-of-magnitude estimates 
for interlopers, we will neglect sub-halos in our analysis. 



2.2.2 The brightest galaxies 

A second way to estimate the number of candidate host 
galaxies is through the M — L relation, where L is the lumi- 
nosity of the host galaxy. Of particular interest is the fact 
that the M — L relation and the M — a relations are dis- 
crepant at the high-mass end (here a denotes the velocity 
dispersion of the host). The former predicts higher masses 
for the most mass ive SMBHs, and h igher number densities 
for fixed BH mass (jLauer et al.ll2007l . and references within). 
This therefore results i n a greater number of individually re- 
solvable PTA sources (jSesana. Vecchio fc Volonterill2009l '). 

Since a is used to infer Afhaio, we expect that for a 
fixed SMBH mass, the number of expected interloping host 
galaxies, inferred from the M — L relation, would also be 
greater than the number of halos, inferred from the M — a 
and a — Afhaio relations. 

To evaluate this different estimate quantitatively, we 



adopt the M — L relation found by (jLauer et al.l 120071 ) for 
the most luminous core galaxies in their sample. 



Mv 



-22.0- 1.8 logioA^g, 



(6) 



where Adv is the V^-band magnitude of the host galaxy. 

To compute the num ber of sufficient l y lum i nous galax- 
ies, we use the results of iGabasch et al.l (|200J, 120061 ). who 
measured the luminosity function in multiple wavelength 
bands between 150 — 900 nm, and studied the redshift evo- 
lution in each band out to 2 > 2. The luminosity function is 
given in the form of a standard Schechter function (jSchechteil 

mi), 



1 ao + l 



c^iMv) = |(lnlO)r [l0('/^)("--«-)]' 

xexp[-10'^/^""^--"-']. (7) 

They set a constant value for the parameter oq while fitting 
My and 0, to a power-law redshift dependence of the form 



Mv{z) = Mv.o + Aln(l + z), 
f(z) = (PU^ + zf. 



(8) 
(9) 



The five fitting parameters (ao, My fj,(j>o, A, B) vary with 
the wavelength ba nd of the luminosity function. Because 
the iGabasch et al.l results do not have fits for the V-band, 
we interpolate the parameters between neighbouring bands 
to obtain the following values: ao 



6.2 X 10 



"^ Mpc"^ 



M-\A: 



=:; -1.3, M,>_o 
-1.18, and B k 



-21.1, 
-1.05. 



2.2.3 The brightest AGN 

Finally, a third method to identify plausible hosts is to 
search for AGN that are luminous enough to be plausibly 
powered by a Af ~ l(f Mq SMBH. AGN activity is an ideal 
scenario for identifying the EM counterparts of PTA sources, 
as the interaction between a compact SMBH binary and its 
accretion flow provide a natural physical mechanism for elic- 
iting a smoking-gun EM signature. However, a significant 
uncertainty with this approach is whether the host of a re- 
solved PTA source is likely to be undergoing an observable 
AGN episode. While multiple studies have suggested that 
galaxy mergers trigger AGN activity (refs. in !JT]) whether 
the two phenomena are causally related remains an open 

question. 

Recently, ISchawinski et al.l (|201ll ) suggested that the 
low Sersic indices in most X-ray-selected AGN hosts at 
1.5 < z < 3 indicate that they are disc gal axies, and there- 
fore unrelated to mergers (see, however, iGovernato et al.l 
I2OO9I , who suggest mergers can result in disc galaxies). Fur- 
ther, even if one accepts that there exists a direct causal con- 
nection between galaxy mergers and luminous AGN activ- 
ity, it is uncertain whether such a trend extends to the most 
massive galaxies at 2; < 1.5. The mass fraction of cold gas in 
massive galaxies tend to decrease toward lower redshift, and 
gas-poor "dry" mergers are thought to play an important (if 
not dominant) role in the assembly of giant ell iptical galax- 
ies at z < 1, in the field as well as in clusters (|van DokkumI 
I2OO5I : iLin etalll2008l . [JOloh . On the other hand, the amount 
of gas required to fuel a luminous AGN episode is a small 
fraction of the total gas content of even very gas-poor galax- 
ies; the most luminous known AGN are situated in giant 
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elliptical galaxies; a plur ality of mergers of massive galax- 
ies at 2; < 1 are gas-rich (jLin et al.ll2008l ): many early-type 
galaxies identified as undergoing a dry merger have been 
found to contain d etectable amounts of gas in foUowup HI 
observations (e.g., Donovan. Hibbard fc van GorkomI 120071 : 



ISanchez-Blazguez et al.ll2009f ): and even though the hosts 
of the most luminous quasars tend to be ellipticals, they 
are not exclusively s o in the SMBH mass regime of inter- 
est (M>10^Mn;e.g., iPercival et al.l200ll : iFlovd et al.l2004 
IZakamska et al.ll200q . and references therein). We conclude 
that PTA sources powering luminous AGN activity is a plau- 
sible scenario, and not merely an expedient assumption. 

We parameterize the minimum luminosity for the 
AGN counterpart in terms of the Eddington luminosity 
L^EddiM) = 47rGM iieTn-pc/oT, where /^e is the mean molec- 
ular weight per electron and ot is the Thomson cross sec- 
tion: 



Lmin(M) = /miniEdd (M) . 



(10) 



We choose /min = 10 ^ for our minimum Eddington ratio 
L/Lfidd, motivated by the fact that t his quantity is observe d 
to peak at L/I/Edd ~ 0.1 - 0.3 (e.g.. lKollmeier et al.lliooel ). 

In order to estimate the number of AGN that are bright 
enough to correspond to a PTA source with mass M, we 
adopt the observationally motivated fits to the AGN lu- 
minos ity function given bv lHopkins. Richards fc HernquistI 
l|2007l : their equations 6, 18 and 20, and Table 3). 

Note that although the Eddington ratio distribution 
and the luminosity function cited above are expressed in 
terms of the bolometric lum inosity, they are actually proxies 
for th e optical luminosity. iHopkins. Richards fc HernquistI 
l|2007h noted that their bolometric luminosity function is ef- 
fectively equivalent to the optical luminosity function, and 
iKoUmeier et al.l (|2006l ) uses the flux at 510 nm to estimate 
the bolometric luminosity. Our AGN interlopers are there- 
fore optically luminous AGN, and we assume nothing a pri- 
ori about the X-ray and UV emission of accreting PTA 
sources. We will discuss the importance of searching for the 
EM counterpart at optical wavelengths in 



2.3 Expected Counts of Interloping Galaxies 

If the individual pulsa r contr ibutions to the signal cannot 
be extracted, as in the lSVlOl scenario, then the chirp mass 
and luminosity distance of the resolved source cannot be 
independently known, and the source can only be localized 
within a solid angle Af2. The only constraints on M and Dl 
are then model-dependent, and come from theoretical expec- 
tations for the population of resolvable sources, given the 
detection threshold of the array. The upper end of the chirp 
mass distribution of SMBH binaries, along with the detec- 
tor sensitivity, sets a maximum luminosity distance Dl, ina.x 
(equivalently, Zmax). Similarly, the chirp mass distribution in 
the local Universe determines a minimum chirp mass Almin 
required for a PTA source to be resolved. Note that since 
M = ri~'-^^^^M > 2^^^M,the quantity Mmin also sets a lower 
limit Mmin ~ 2.3A1min ou the gravitational mass. 

The number of interloping halos can be expressed as 



iV, 



(SV) 



An 

4tt 



A/halo(A/„ 



dnhaio ,,, dy , 

-j—f dMhaio -J- az, 

dAfhaio az 

(11) 




Figure 1. Estimates of the number of interloping host objects 
— (a) massive dark matter halos, (b) luminous galaxies, and (c) 
luminous AGN — in the conical error volumes suggested bv lSVlOl . 
The extent of the error volume is limited by zimax, the maximum 
redshift at which PTAs can resolve an individual source, and the 
angular localization AC = 40deg . The number of interlopers 
is calculated by assuming a minimum SMBH mass Mmin, which 
then sets the minimum host mass/luminosity through equations 
\E\M andlTol 



where n is the comoving number density of dark matter 
halos, and dV/dz = 4,-kD\ dDh/dz is the comoving volume 
element. The lower limit of the integral over halo mass is 
given by equation [5] Similarly, the numbers of interloping 
galaxies and AGN are given by the expressions 



N. 



(SV) 



iV 



(SV) 
AGN 



Ml 

'^ JO 

An 

47r 



dMv -T- dz, (12) 
-oo dz 

^dL^dz.(13) 

'min(A-fmin) 



The limits of integration over luminosity are taken from 
equations [6] and [10] 

We show in Figure [T] the estimated n umber of inter- 
lopers for the worst-case error box in the ISVld scenario, 
assuming An = 40 deg'^ , as a function of the maximum red- 
shift Zmax and minimum BH binary mass Mmin. Panels (a), 
(b) and (c) show the isonumber contours of the expected 
number of interloping massive halos, luminous galaxies and 
luminous AGN, respectively. All three methods to estimate 
the number of interlopers yield on the order of Ng > lO'^ 
for PTA sources with M > 1O'A/0, if the redshift range is 
restricted to Zmax ~ 1. 

The difference in the number of interlopers between the 
top two panels (halos vs. galaxies) for the most massive 
SMBHs arises because the observed SMBH samples yield 
an internally inconsistent set of M — cr, L — a and M — L 
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Figure 2. Same as Figure [T] exc ept that the error volume is 



calculated from the results of lCClCl who assumed that the pulsar 
term of the GW signal can be used to infer the luminosity distance 
to the source binary. The error box is limited by the uncertainty 
ADi^/Di^ = 20% in the luminosity distance to the source, and 
the angular localization Af2 = Sdeg'^. Note that whereas the 
horizontal axis in Figure [T] showed the maximal PTA detection 
range ^max, here it denotes the actual redshift z of the source. 



relations, as mentioned above. While the interpretation of 
this inconsistency is bey ond the scope of our paper, we note 
that lTundo et al.l (|2007l ') discussed this issue, and concluded 
that the intrinsic scatter in the relations produces a selec- 
tion bias: using the observed BH samples yields a biased 
L — a relation (too low L for given a). This suggests that 
the M — L galaxy relation we adopted may also be biased 
and it under-predicts L; correcting this bias would decrease 
the number of galaxy interlopers. 

If the GW signal can be used to constrain A4 a nd Dl 
of the source via statistical inference, as suggested bv lCClOl . 
then the numbers of interloping halos, luminous galaxies and 
AGN are given by 



An 

4-K 






duh 



dMh 



dMhalo,(14) 



Mhalo(M„i„,z) 

AQ. f~~+ f°° driAGN 

47r 



(/> dMv ^ dz, (15) 
dz 



imi„(J^-fmi„) 



dL 



dL, 



(16) 



respectively. Above, the redshifts z± = z(D_l± ADl) bound 
the radial extent of the error box. We adopt AQ = 3 deg^ 
and ADl/Dl = 20%. We ignore errors due to weak 
lensing, which are expected to be on the order of sev- 
eral percent for sources with z ^ 1.5 (JKocsis et al.l 120061 : 
iHirata. Holz fc Cutlei|[201ol : IShang fc Haimanll201ll ) We do 
not place an upper limit on the host halo mass (or on the 



host galaxy luminosity). In principle, such an upper limit 
could be computed, given PTA's observational error on the 
chirp mass and the spread in the ratio between the chirp 
mass and the gravitational mass of the binary (i.e., from 
the model-dependent mass ratio distribution of resolved 
sources). For example, [CCly provide a chirp mass error esti- 
mate of AA1 ~ 5%. Converting the chirp mass to the grav- 
itational mass, however, can introduce a large uncertainty, 
e.g. a factor of ~ 2 depending on whether the mass ratio 
is 0.1 or 1. Since the number density of interlopers decrease 
rapidly with increasing halo mass (or luminosity), this sim- 
plification should not affect our estimates. 

In Figure O we plot the number of interloping host 
candidates against the source redshift z. Not surprisingly, 
the prospe cts for EM identification improve dramatically in 



prospe 

Iccid 



the ICClOl scenario. For massive (M > 10 Mq) resolved 
PTA sources, we anticipate that the error box will con- 
tain a single host candidate at z ^ 0.2, and several hun- 
dred at z ^ 0.7. We expect only a single group-sized halo 
(M > few X 10^"^ M0) in the error box at any redshift in the 
ICCld scenario. Note that the number Ng of interlopers is 
not necessarily a monotonically increasing function of 2, as 
the decline in the number densities of the interloping objects 
(in particular massive halos) competes with the increase in 
the comoving size of the error boxes. 

Our simple calculations show that in the scenario of 

ICCld . resolved PTA sources with M > lO^Mg and z <, 0.5 
are likely to have at worst dozens of interlopers in the error 
box. With this low number, one could conceivably perform 
follow-up observations of each individual candidate. If, on 
the other hand, luminosity distances to the source cannot 
be determined, this number increases to ~ 10^, suggesting 
that it will become extremely difficult to electromagnetically 
identify the source in the absence of an obvious, tell-tale EM 
signature. 

In practice, the number of interloping galaxies may 
be somewhat larger than the value computed by equations 
I11H 16I The halo mass of any given candidate host system 
will not be known a priori, and the intrinsic scatter in the 
AfsMBH — Mhaio (MsMBH — o-hoat) relation will lower the min- 
imum halo mass threshold for candidacy. On the other hand, 
the simple calculations presented here do not consider de- 
tailed demographic properties of resolved PTA sources and 
plausible host s, such as the presence of a nucl ear stellar 
core (iMakino Il997: iRavindranath. Ho fc Filippenko 2002 



iMilosavlievic et al.ll2002l : IVolonteri. Madau fc Haardtll2003[ ) 



or galaxy morphology. Including such factors in the analysis 
will narrow the field of candidate hosts. 



As we argue in 13.31 candidate AGN counterparts may 
be further vetted by examining their UV and X-ray emis- 
sion fo^Jeatures indicative of a central SMBH binary (see 
also lSesana et al.ll201ll for an in-depth discussion of possible 
high-energy signatures for pre-decoupling — i.e., tew > t^, 
— PTA sources). In addition, PTA sources are sufficiently 
nearby that it should be possible to observe an interloping 
AGN together with its host galaxy. It should therefore be 
possible to combine the AGN emission, the galaxy luminos- 
ity and the inferred SMBH mass to cross-check candidate 
counterparts. 
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3 ACCRETION DISCS AROUND 
PTA-SOURCE BINARIES 

Motivated by the results of the previous section that the 
number of plausible host galaxies in the PTA error box 
may be tractable for follow-up EM searches, we next model 
the EM emission properties of SMBH binaries detectable 
by PTAs. We focus our attention on SMBH binaries that 
are undergoing luminous accretion, as these are the most 
promising class of objects for EM identification. 

Normalizing the binary mass M and rest-frame period 
P to the typical orders of magnitude expected of resolved 
PTA sources, M = lO^A'/oMg and P = 1 yr Pi, we write 
the semi-major axis for the source binary as 



a{M,P) 



-2/3 p2/3 



GM 



2.23 X IQ-'^Mg^^P^^''^ pc 



(17) 



Binaries detectable by PTAs have long overcome the so- 
called "final parsec" problem. The rest-frame time to merger 
for a binary with mass M and semi-major axis a, driven by 
GW emission alone, is 



256 G^M^ 



V 



1.9 X 10^ Mgr?!"! 



yr 



(18) 



102GAf/c2 
= 2.0 X 10^Mg"^''^77j:lPf''^ yr 

(jPetersI 1 19641 ). Because typical resolved sources have 
M/M = r]~''''^ ~ 3, we normalize the symmetric mass 
ratio T] = {M2/Mi)/[1 + M2/Mi]^ to the value r;i:4 = 
?7(M2/Mi = 0.25) = 0.16. Note that our ad hoc transla- 
tion between M and M is not very sensitive to the value 
of q; the ratio M/M varies by less than a factor of two in 
the range 0.1 < M2/M1 < 1. The upper bound in equation 
[18] corresponds to binaries in circular orbits, with eccentric 
orbits merging faster. Recent work has shown that bina- 
ries may have eccentricities a s high as ~ 0.6 at decouplini 



(IRoedig et al 



ICuadra et al. 



20 111 : see also lArmitage fc NataraiarJ |200_. 
II2OO9I ). Thus, typical PTA-resolved sources will 
coalesce on scales of ~ 10'^ years. However, exceptionally 
compact sources will coalesce on scales of several years; for 
example, a binary with P = 0.1 yr — approximately the 
lowest binary period that is expected to be observable with 



PTAs — will merge in fn 



' 4 yr. 



The tidal torques of the compact SMBH binary 
provide a particularly promising mechanism for producing 
a tell-tale observable feature. Theoretical calculations 
jGoldreich 



Tremaine 



198d: Artvmowicz et al.l 



Artvmowicz fc Lubowl 1994: lArm itagc & Nataraian' '200 
Bate et al.l 12 003: Havasaki. M incshiec & Sudou 20£ 
MacFadven fc Milosavlievid I2OO8I : ICuadra et al.l I2OO 



Chang et al.l I2OI0I ) robustly predict that in geometrically 
thin circumbinary accretion discs, binary torques can 
open an annular, low-density gap around the orbit of 
the secondary. The gas inside the gap accretes onto the 
individual SMBHs while the gas outside is pushed outward 
by the tidal torques. The binary's tidal torques transfer 
orbital angular momentum into the outer disc, causing 
the binary's orbit to shrink gradually while maintaining a 
roughly axisymmetric circumbinary gap. 



T he gap opens near the re sonance radius R ~ 3^ "a ~ 
2.08a (|Artvmowicz et al.lll99ll ). Numerical simulations of 
thin circumbinary discs (see refs. in above paragraph) pro- 
duce gaps with an azimuthally averaged radius of 1.5 — 3 
times the binary semimajor axis. The exact size and shape 
of the gap is not easily characterized; the geometry depends 
on the binary masses and orbital eccentricity, as well as the 
efficiency of angular momentum tr a nspor t within the disc. 
Following iMilosavlievic fc PhinnevI (|2005|), we parametrize 
the size of the gap as R\ = 2Xa, where A ~ 1 is a di- 
mensionless parameter. We are interested in circumbinary 



discs that are truncated inside _Rx 



200Mc 



-2/3 p2/3 



P^'-'GM/c-' 

(equation 1 17|) . Below, we model surface density profiles and 
thermal emission spectra of such discs, as well as the ther- 
mal emission due to leakage of gas into the cavity and onto 
individual SMBHs. 



3.1 Disk Properties and Binary Decay 

3.1.1 Disk around a solitary SMBH 

Adopting a geometrically thin, thermal gray-body disc 
model (e.g., lBlaesll2004 IMilosavlievic fc Phinnevll2005h , we 
estimate the properties of circumbinary discs around re- 
solved PTA sources. As a reference model, let us consider a 
disc around a solitary SMBH. 

Until recently, the standard a-viscosity prescription 
for accretion discs, in which the kinematic viscosity scales 
proportionally with the total pressure in the fluid (i.e., 
f oc a(pgas + Prad)), was thought to be thermally 
and viscously unstable (|Shakura fc Sunvaevlll976l : IPringlel 
1 19761 ). Altho ugh magnetohydrodynamic s hearing box sim- 
ulations by iHirose. Krolik fc BlaesI |2009| have since sug- 
gested that a-discs are actua lly thermally stable, they 
may still be viscously unsta ble (|Lightman fc Eardlevlll974l : 
iHirose. Blaes fc Krolikll2009r i. We therefore assume a viscos- 
ity prescription in which the kinematic viscosity scales with 
gas pressure (a.k.a. the "/3-disc" model): 



_ 2 apgas _ 2 aksT 
3 pO, 3 firripQ 



(19) 



This viscosity prescription is consi stent with previous anal- 
yses of thin circumbinary d iscs (|Milosavlievic fc PhinnevI 
I2OO5I : iTanaka fc Menoull201Ct l. 

If discs around PTA sources are instead described by 
the standard a viscosity prescription, then they would have 
somewhat higher temperatures than what we calculate be- 
low using the /3 prescription. The higher viscosity in the a 
model would also result in the gas being able to follow the 
decaying binary to closer separations. Both effects would re- 
sult in the thermal spectrum being somewhat harder than 
we calculate below, but we expect our qualitative results to 
hold as long as the disc is able to remain geometrically thin. 

Observations of accret ing Galactic compact objects 
King. Pringle fc Livid 120071 , and refs. therein) and blazars 



Xie et al. 



2009h 



suggest a ~ 0.1 — 0.4, while studies of 
AGN continuum var iability place a lower limit of «>0.01 
jStarling et al.l 120041 '). Numerical simulations are consistent 
with a > 0.1, with lower reported v alues possibly being due 
to sm all sizes of the simulation box (|Pessah, Chan fc PsaltisI 
[20071). We choose a = 0.3 as the fiducial value, and write 
00.3 = a/0.3. We assume that shear viscosity is the domi- 
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nant mechanism of radial gas transport in the circumbinary 
disk. 

The surface density E and the mid-plane temperature 
T of the disc are obtained through the following equations: 



l{Q.,Tp)as^T^ — 



p 

'p 

r 



8 
3r 

M 



(20) 

(21) 
(22) 

(23) 



Above, H is the deviation of the bolometric flux from black- 
body due to the photons be ing thermalized above the mid- 
plane (see, e.g., [Blaeg[200J), Tp is the temperature of the 
thermalization photosphere, and r is the optical depth be- 
tween the mid-plane and the thermalization photosphere. 
The quantity ^ is a porosity factor that relates the surface 
density to the optical depth. We set it to 0.2 following [Turner] 
l)2004h and express our results in terms of 9o.2 = 9/0.2; 
however, most of the disc properties are not very sensitive 
to this parameter. We will use the dimensionless parameter 
rh = M /M-Edd to describe the local accretion rate in units 
of the Eddington rate, assuming a radiative efficiency of 0.1, 
i.e. LEdd = O.lMEddC^. 

The disc scale height H is evaluated in the usual way: 



P 



(24) 



where Cs = \/p/ p is the isothermal sound speed calculated 
from the total pressure p = pgas -I- Prad, and the volume 
density p of the disc is given by p = T./ H. The second term 
on the right-han d side of equati on 1241 is due to the disc's 
self-gravity (e.g.. iPaczvuskiHigTSl ). 

In the regions of interest, the primary source of opac- 
ity is electron scattering, and radiation pressure dominates 
over gas pressure. In this regime, the gray-body factor S 
can be approx imated as H Ri 0.17(n yr)^/^[Tp/(10*K)]-^^/^ 
(jTanaka fc Mc nou 2010|), and with a little algebra we obtain 
the surface density profile in the disc: 



E(i?) ^ 1.3 X lO'^gcm^^ ' 



R 



xM, 



V lOOGM/c 

,16/85 .36/85 -4/5^-1/5 



^0.3 ''0.2 



-6/17 



(25) 



A steady-state disc far from the central object satisfies M = 



SnuT, — constant, and so we have u ocT, 



E"/! 



3.1.2 Circumbinary discs around orbit- decaying binaries 

After a circumbinary gap is opened, the SMBH binary un- 
dergoes several stages of orbital decay. Let us briefly examine 
the different stages, and the orbital evolution timescale (or 
residence time) frea = a/\da/dt\ for each. Our goal here is to 
describe the structure of a dense gaseous annulus, extending 
at least a factor of few in radius, that is created around the 
PTA source. The annulus results from inward migration of 
the binary from larger radii in a more extended accretion 
disc. For a more thorough discussion of the orbital decay of 
SMBHB bina ries, through various physical reg imes in a thin 
disc, see, e.g. iHaiman. Kocsis fc Menoul ((20091). 

We begin with disc-driven orbital decay, in which the 



binary's tidal torques transfer its orbital angular momen- 
tum to the surrounding gas. At large orbital separations, 
the mass of the gas at the edge of the cavity far exceeds 
the mass of the secondary. In this regime, analogous to disc- 
dominated Type II migration for proto-planets, the binary's 
orbital evolution is limited only by the rate at which the 
nearby gas can transport away angular momentum, i.e. 



.(disc) 



tu{R\ 



2Ri 



(26) 



The tidal torques prevent the gas from flowing inward of 
Rx, and so the region inside the gap is starved. Any gas 
that is ini tially present will b e depleted on the local viscous 
timescale (jChang et al.l20ld ). In standard steady-state thin- 
disc models the viscosity is an increasing function of radius, 
so this drainage occurs on timescales shorter than that of 
the binary's orbital decay. 

When the mass of the secondary becomes comparable 
to the local disc mass, the orbital decay slows down with 
respect to the local viscous time. The gas piles up imme- 
diately outside the cavity, forming a decretion region in 
which the vi scous torque 71- — 2i-kvT,Q,R? is nearly constant 
with radius (IPringle lOQlll. We apply the analytic model of 
llvanov. Papaloizou fc PolnarevI (|l999l ) to calculate the resi- 
dency time for this secondary-dominated migration stage: 

,(scc) _ yM 



-tiy{Rx 



(27) 



inRlE{Rx) 

Note that there are two competing effects influencing ti-H'^' : 
the decay slows down as the local disc mass decreases 
with respect to the secondary, but this is mitigated to 
a small extent by the fact that E outside the cavity in- 
creases due to pile-up. The enhancement of E relative 
to that of a disc around a solitary SMBH of the same 
mass as the binary (equation | 25[l ha s the functional form 
(jlvanov. Papaloizou fc Polnarevll 19991 ) 



j-j(binary) 
J](aolitary) 



= <^ 1 + A 



1- 



Rx 



R 



(disc/sec) 
A 



1/2 



R_ 
Rx 



^1/2 



in the neighbourhood R > Rx. Above, R\^ isc/sccj 



(28) 
is the ra- 



dius of the cavity when the transition from disc-dominated 
to secondary-dominated migration occurs, i.e. when ryAf = 
A'jTRf_'S{Rx). For reasonable parameter values, R^^ isc/aoc) ^ 
10 GM/c . The dimensionless quantities A and B in equa- 
tion [28] depend on the viscosity and mass p rofiles of the 
disc (see llvanov. Papaloizou fc PolnarevI Il999l . for details) . 
We typically find that A ~ 4 and B ~ 0.2 in our disc 
models; i.e., the fractional surface density enhancement 
during secondary-dominated migration is no greater than 
(1 + A)^~ 1.4. 

At yet smaller separations, the binary's orbital evolu- 
tion begins to be driven by GW emission. Since binaries 
of interest here are far from merging, GW emission can 
be approximated by the leading term in the Newtonian 
qu adrupo l e. For circular orbits, the residence time is given 
bv lPetersI (|l964l ) 



,(GW) _ . 

•res — ^(rncrgc 



Jl_l^ 1. ('291 

64 G3M3 r? ■ ^ ' 

As the binary's orbital decay accelerates due to GW emis- 
sion, the pileup caused by secondary-dominated migration 
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spreads out. Past the point where frca ~ tv{R\), the bi- 
nary begins to outrun the disc, as the decay timescale for 
a becomes rapidly shorter than that on which the disc can 
viscously spread. 

Let us now discuss the gravitational stability of the disc, 
based on the stability crite ria of a radiation-p ressure domi- 
nated fluid summarized by iThompsonI (|200a ) . If the radia- 
tive diffusion timescale is much shorter than the dynamical 
timescale, then the radiation pressure does not stabilize the 
fluid and gravitational fragmentation occurs on the same 
length scales as it would in the absence of radiation pres- 
sure. If the radiative diffusion timescale is much longer than 
the dynamical time, which we find to be the case for our disc 
models, then radiation pressure acts to make the fluid more 
Jeans-stableQ We invoke the Toomre criterion, and assume 
that the disc is gravitationally stable when the dimensionless 
parameter 



Q{R) 



ttGE 



(30) 



is greater than unity. Note that the key effect of radiation 
pressure in this context is that the sound speed Cs is com- 
puted from the total pressure, not just the gas pressure. 

A counterintuitive result is that in radiation pressure- 
dominated discs, increasing the surface density — or equiva- 
lently, the accretion rate — at fixed radius and disc parame- 
ters (a, 6, etc.) makes it more stable against fragmentation. 
This is a significant point, because it allows for the exis- 
tence of a copious amount of hot, gravitationally stable gas 
near the binary. We demonstrate this behavior as follows. 
Combining equations [24] and 1301 we may write in general 

Q 



- W: 



x^ — X, where x = p/{4ttGT, ) > 1 is the ratio of the 
pressure in the disk to its self-gravity. It follows directly that 
the stability criterion Q > 1 is equivalent to the condition 
a; > (2 -I- %/5)/4 ~ f.f . Some algebraic manipulation of the 
disc equations (|20l — I24p gives another general expression 
for the region dominated by electron scattering and radia- 
tion pressure: 



X — —il v^ oc — . 



8.CG E- (31) 

In standard a and /3 discs, increasing the surface density 
with all other properties held constant increases the mid- 
plane temperature. This raises the value of y, and the gray- 
body factor E decreases or remains constant (e.g., H = 1 
for blackbody discs). Thus, in radiation pressure-dominated 
discs the parameters x and Q increase with an increase in 
E or the accretion parameter m. Increasing the accretion 
rate heats and "puffs up" the radiation pressure-dominated 
portion of the disc more than it acts to enhance self-gravity. 
The opposite is true in the gas pressure-dominated region, 
where one obtains x ex 'il~^'^H^^K^'^'S.^^'^\ increasing the 

The gas is susceptible to an additional weak diffusive instability 
that grows on the Kelvin-Helmholtz timescale, Jkh ~ K,c^/{nGc). 
We find that the viscous timescale is shorter than the Kelvin- 
Helmholtz timescale — i.e., the diffusive instability is irrelevant 
— in all but the outermost annulus of the radiation-dominated 
region of our fiducial circumbinary discs (e.g., at R = 400GA'//c^, 
^KH ~ 2 X 10*^ yr and f^ ~ 10^ yr). Even in the small region 
where the disc is formally unstable to the diffusive instability, 
it is plausible that local turbulence can quench its growth. We 
therefore assume the diffusive instability is unimportant. 



surface density (or accretion rate) in the outer disc drives it 
closer to gravitational instability. 

We find that the disc is Jeans-stable inside a radius 



Rq 



R{Q = f ) 






•'o.s 



Mc 



-16/55fll7/55 
''0.2 • 



(32) 



Equation [32] is valid for radiation pressure-dominated re- 
gions only. As a function of the accretion rate m, Rq has 
a minimum value, where it coincides with the radius where 

Prad — Pgas • 

The gas density profile in the outer regions R > Rq, 
where classical thin-disc models predict Q < f , is uncertain 



One p ossibility that has been explored bvlSirko fc GoodmanI 
12003") and others I^Thomp son. Quataert fc Murravl l2005l : 
[Levin 20071: iLodato et al. 2003) is that feedback mechanisms 
(such as nuclear fusion from stars that formed in the disc 
or their supernovae) inject sufficient energy as to maintain 
marginal gravitational stability with Q ~ f in the outer re- 
gions. However, the profile of the outer disc is not central to 
this study, as we are interested in radiation from the central 
regions of the disc, where the presence of a compact binary 
is most likely to produce characteristic features that may 
distinguish them from accretion discs around single SMBHs. 
To keep our analysis as simple as possible, we simply neglect 
the thermal radiation of the disc outside Rq . 

What is the accretion rate in the disc? Uncertainty re- 
garding the outer gas distribution notwithstanding, quasars 
are able to efficiently supply SMBHs of mass > 10* M© with 
enough fuel to maintain luminosities of 0.1 — Il/Edd for pe- 
riods of 10®~* yr. That most quasars radiate at just un- 
der the Eddingt on limit while few exceed LEdd (see, e.g., 
[Shen et al.|[2008l . for quasar Eddington ratios at z ^ 1) sug- 
gests that their luminosities are limited by radiative feed- 
back, rather than by the availability of fuel. If this is the 
case, then the surface density in a discs around compact bi- 
naries can be significantly greater than in a disc around a 
solitary SMBH of the same total mass. This is not because 
of the mass accumulation of gas outside the binary's orbit, 
but because such a disc would have a much lower luminosity 
owing to its low-density central cavity. If binary torques in- 
hibit gas from accessing the centre of the potential, then this 
effectively reduces the radiative efficiency of the system, i.e. 
idisc ^ Mc^. We find that even at m ~ 10, the locally vis- 
cously dissipated flux in our discs does not provide sufficient 
radiation pressure to unbind gas from the local gravitational 
fieldQ 

In Figure [3] we plot, for a circumbinary disc around 
a binary with Afg — 1, M2/M1 = 1/4, several transition 
radii as a function of the mass accretion parameter rh. We 
show the radii at which the disc transitions from being radi- 
ation pressure-dominated to gas pressure-dominated; from 
where the opacity is dominated by electron scattering to 
free-free absorption; and from Jeans-stable to unstable. We 
note that in general, radiation pressure acts to stabilize the 
disc against Jeans collapse, and that the radius Rq closely 



"^ We remind the reader that rh is defined with respect to the 
Eddington limit assuming a radiative efficiency of ~ 0.1. Strictly 
speaking, our circumbinary discs are not super-Eddington, oven 
when the parameter rh exceeds unity. 
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Figure 3. In black lines, several transition boundaries within 
a steady-state, thin accretion disc are plotted as a function of 
the mass supply rate rh = M/{0.1L^^^/c^). Binary parameters 
are M = IO^Mq, M2/M2 = 1/4; disc parameters are a = 0.3, 
X = 1, 9 = 0.2. Plotted are the radii in the disc where the 
disc is marginally stable to gravitational fragmentation [Rq; dot- 
ted lines); where the radiation pressure equals the gas pressure 
(dashed); and where free- free opacity equals electron-scattering 
opacity (dash-dot). In blue lines, we plot the size of the cir- 
cumbinary cavity Rx = 2Xa where the binary's orbital decay 
transitions from being gas-driven to GW-driven (thick lines), and 
when the orbital decay timescale becomes shorter than the viscous 
timescale at the inner edge of the disc (thin lines). For reference, 
the radius of the cavity for a 10® Mq binary with a rest-frame 
period of 1 yr is plotted as a horizontal red dashed line. 



corresponds to the radius where the disc becomes radia- 
tion pressure-dominated. We also plot the size of the cavity 
R\ ~ 2a at which the binary's orbital evolution transitions 
from being gas-driven to GW-driven, and the value of R\ 
where tf^ = t^{R\). We find that the disc is geometrically 
thin {H/R < 1) or marginally thin {H/R ^ 1) for m < 10. 
In conclusion, Figure [3] suggests that a Jeans-stable cir- 
cumbinary annulus could exist, instantaneously, around an 
individually resolved PTA source, for any value of the sup- 
ply rate, extending at least by a factor of two in radius (from 
the inner radius shown by the [red] dashed line, to the outer 
radius shown by the [black] dotted curve). However, in order 
for this annulus to be created through the in-ward migration 
of the secondary BH from larger radii in a more extended 
disc, we require that the disc is stable to radii that extend 
beyond the gas/GW-driven transition. This latter require- 
ment (i.e., the dotted [black] curve must lie above the thick 
[blue] solid curve) means that in practice, the gaseous an- 
nulus exists in PTA sources only if the mass supply rate is 
TO > 1. As argued above, while radiative feedback may disal- 
low such high (super-Eddington) rates in thin discs around a 
single BH, they can naturally be maintained in discs around 
binaries, owing to their low radiative efficiency. 



3.2 Surface Density Evolution of the 
Circumbinary Gas 

Let us now address the surface density profile of the outer 
disc at the time when the binary SMBH becomes observable 
by PTAs. 

The tidal torque density dTtidc/dR is sharply peaked 
in a narrow region that roughly coincides with the edge of 
the cavity R\, preventing the gas from accreting inward. 
Everywhere else in the disc, the tidal torques are negligi- 
ble compared to the viscous torques. The effect of the tidal 
torques in the region R « Rx can thus be approximated 
as a boundary condition prohibiting mass flow ac ross R\ 
(|Pringldll9"9ll : llvanov. Papaloizou fc Polnarevlll999l '): 



M{Rx,t) = Qtvi^T, 



dlnjEuR 
dlnR 



1/2N 



0, (33) 



Note that our disc is not steady-state, and the local mass 
flow rate M need not be radially constant. 

The surface density evolution of the circumbinary disc 
Rx is governed by the standard equation for viscous discs 
(e.g.. |Pringi3ll98ll : [Rank. King fc Rainell2002l ) without in- 
cluding an explicit term for the tidal torques: 



i^f''-'^--nm 



-,1/2 



d 



i.--(3.Ei? 



jl/2 



(34) 



A semi-analytic solution for the thin-disc equation 1341 
with th e boun dary condition in equation [321 was derived by 
iTanakal ((20111), for a flnite boundary Rx > and a special 
viscosity prescription v oc R". The solution can be written 
in the form 



/"OO 

E(i?, t)= G{R, R'; t, Rx) Einit(i?') dR' , 

J Rv 



(35) 



where G{R,R';t,Rx) is the Green's function specific to the 
boundary condition and the chosen value of the viscosity 
power-law index n; and Einit(-R) is an arbitrary initial den- 
sity profile. We find that our discs satisfy v oc 7?°* inside 
R < lO^GM/c^, and thus adopt n = 0.4. 

In order to model a thin accretion disc around a GW- 
driven SMBH binary, we modify the solution of iTanakal 
(|201H ) in two ways. First, we derive a more general Green's 
function to allow for a boundary condition with nonzero 
mass flux across the inner boundary: 



M(iiA,f) = /loakM,,(i?A,f). 



(36) 



Above, Mas = SirvT, is the standard accretion rate for ex- 
pected of a steady-state disc, and < /icak < 1 is a numer- 
ical factor representing the incomplete suppression of gas 
inflow into the cavity. The case /icak = corresponds to 
total suppression of accretion by the binary's tidal torques. 
The boundary condition in equation [31] is motivated 
by results from numerical simulations, which show that 
in general the binary torques do not completely pre- 
vent accretion into the gap, but rather allow some gas 
to leak into the centre of the disc with a suppressed 
mass flux /leak ^ 0.1 (e.g., Artvmowicz fc Lubow 1996|; 

lOchi. Sugim oto fc Hanawal 



Giinther. Schafer fc Kle^ 12004 



20051 : iMacFadven fc MilosavlievidbOOSl '). We choose /icak = 



0.1 as our fiducial value. 

The long-term behavior of the gas is to pile up near 
the cavity and satisfy the power-law i/E oc /^(Jitjak-iJ/a -^^ 
the vicinity of the boundary. In comparison, a steady-state 
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around a solitary central mass disc satisfies i^E = constant, 
and a boundary condition imposing zero inward mass flux 
satisfies i/E oc _R~^' at the boundary. Note, however, that 
the value of /icak does not have a strong effect on the 
mass profile (and hence the luminosity produced) outside 
the cavity-opening radius Rx ; the fractional surface density 
enhancement due to secondary-dominated migration is typ- 
ically of order unity. 

T he gas that enters the cavity does so in nearly radial 
orbits iMacFadven fc Milosavlievid ()2008l ). and so is dynam- 
ically decoupled from the circumbinary disc. Thus, the sur- 
face density and mass flux inside R\ can consistently be 
disregarded in the Green's function formalism. The leaked 
gas can p resumably form accretion discs around one or both 
SMBHs (|Havasaki Mineshige fc Sudoul 120071 ') . presumably 
at the usual AGN radiative efflciency ~ 0.1. The mass sup- 
ply rate of such circum-secondary (or -primary) discs will 
be modulated by /icak^fss, which decreases as the binary 
outruns the circumbinary gas. Because the viscous time at 
the outer edge of such discs are shorter than at R\, they 
will be nearly steady-state, with the surface density profile 
at any given time being determined by the instantaneous 
mass flux into the cavity. Thus, the bolometric luminosities 
of the discs around each disc may be roughly expressed as 
L < 0.1 /icakMss(-RA)c^, and would be Eddington- limited 
by the potential of the individual black holes they orbit. 
This suggests that if the quantity fica.k'rn exceeds unity, the 
region inside the cavity would develop radiation-driven out- 
flow winds. Thus, the parameter /leak affects the energetic 
output due to accretion inside the cavity far more than it 
does that of the circumbinary disc. 

The Green's function for the boundary condition in 
equation [3n] is given by 



o(-rt, R ; f , Rx 



= '1-2'^" 



-1/4^/5/4 



[Mky)Y,ikyx) ~Y,{ky)Mkyx: 
X \^Mky')Ye{kyx)~Ye{ky')Mkyx: 
X [j!{kyx) + Yt\kyxf 



X exp [-3Afc tl k dk. 



(37) 



Above, £ = 1/(4 - 2n) and A = (1 - n/2fuR-" are 
constants. We have introduced the variables y = R^^"^^, 



y' = R'^ "' ^ and yx = R^^ " , as well as the functions 



Mx) 



X Ji-i{x) ■ 



2-n 



Ji{x) 



and 



xY,.^{x)-^Y{x). 



(38) 
(39) 



2-n 

Taking /icak -^ leads to the solution given in lTanakal (|201ll : 
his equation 42) imposing M{R\) — 0. 

The second modification to the Green's function for- 
malism is to allow the inner boundary to move inward as 
an explicitly known function of time, i.e., Rx{t) = 2\a{t). 
This is done through a time- weighted superposition of differ- 
ent Green's functions at intermediate values of Rx{t). The 
"master" Green's function Q for a time-dependent boundary 
takes on the form 



giR,R'-t) = / dG{R,Ry,Rx) dR>^ ^^ 



dRx 



dt 



(40) 



A derivation and expression for Q is provided in the Ap- 
pendix. 

We take our initial condition as the disc profile when 
^ros^ ~ tv{R\), just as the binary is just beginning to out- 
run the circumbinary gas. We note that the torques exerted 
by the disc are not entirely negligible at this stage. For sim- 
plicity, we approximate the contribution of the disc torques 
to the orbital decay rate, a/{da/dt)d\sc = ircs'^ , as being con- 
stant. This is justified as follows. The disc torques are a weak 
function of _Ra, at least as long as the quantity v{Rx)'i^{Rx) 
is comparable to the steady-state value. Once GW emis- 
sion dominates the orbital decay, the contribution of disc 
torques becomes quickly negligible regardless of the value of 
u{Rx)T.{Rx). 

The orbital decay is then given by 



da 
IE 



da 

GW \ ™ / disc 



a 64 C 

(sec) ^ 



n 



5 G^M^ a3 



(41) 



With the assumption that ires is roughly constant, this has 
the analytic solution 



.(sec) 



t 



In 



a;5-h64c-"^r?tS"'/(5G3M3) 



(sec) 



a4(i) + 64c5r?Cr7(5G3M3) 



(42) 



In Figure |4l we plot a model surface density profile of 
the circumbinary disc around a binary with Mg = 1 and 
M2/M1 = 1/4. As the initial condition, we take a steady- 
state surface density profile with rfi = 3, at the time when 
i{GW) = iv{R\) (solid black curve) |j Initially, the disc has 
a cavity radius of Rx ~ blQGM/c? and a rest-frame period 
of P = 11 yr. We then evolve the profile using our Green's 
function to when the orbital period is P = 1 yr (short- 
dashed blue curve) and P = 0.1 yr (long-dashed red curve). 
We denote with a green dot-dashed line the radius where 
(5=1, beyond which the disc is expected to be susceptible 
to the Jeans instability. All the disc profiles are truncated 
at twice the binary's semi-major axis, and lose mass across 
this radius through the boundary condition in equation 1361 
with /leak = 0.1. Note that the boundary radius Rx moves 
inward faster than the gas can pile up. We see that a small 
amount of gas is able to follow the binary's orbital decay, 
even though the bulk of the circumbinary disc is getting left 
behind by the inspiraling binary. 



3.3 Thermal Emission of Accreting PTA Sources 

Since the innermost gas is missing from the accretion discs 
around PTA binaries, it is probable that their accretion fiows 
will emit less UV and thermal X-rays compared to ordinary 
AGN powered by solitary BHs. The question then is how 



^ Strictly speaking, the surface density profile at this time should 
deviate somewhat from the steady-state one. During secondary- 
dominated migration, the circumbinary surface density profile 
T,{R>^Rx) can become greater than the steady-state profile by 
at most a factor of ~ 1.4 (equation (28] note that the pileup must 
be smaller if one accounts for the fact that /leak > 0). Prior to 
our initial condition, GW emission accelerates the binary's orbital 
evolution, and the circumbinary pileup spreads out; however, the 
surface density does not decrease below the steady-state profile, 
since tres > tv{Rx)- A difference in S of less than 40% is in- 
significant compared to the other theoretical uncertainties, and 
we employ the steady-state profile for simplicity. 
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R (GM/c^!) 

Figure 4. The surface density profiles S for a circuinbinary disc 
around a PTA source. The binary's mass is 10® Mq and its mass 
ratio is M2/M1 = 1/4. We adopt a moderately high value of the 
accretion parameter m = 3, and assume that the circumbinary 
gas can leak into the cavity at the rate given in equation 1361 
with /leak = 0.1 . The solid black curve shows the surface den- 
sity profile when GW emission begins to dominate the binary's 
orbital decay (U(Rx) = taw, Rx = 510GM/c^, P = 11 yr). 
Using the semi-analytic method described in the text, we solve 
for the surface density profile in the disc at later times, when 
P = 1 yr (short-dashed blue line) and P = 0.1 yr (long-dashed 
red line). The dot-dash green line denotes the radius inside which 
the circumbinary disc is stable against Jeans collapse. 



UV- and X-ray-deficient these objects are; the answer de- 
pends on how much gas is able to follow the binary's decay- 
ing orbit, and how much of this gas is further able to leak 
into the cavity and accrete onto individual SMBHs. This is 
a complex problem characterized by dynamical richness in 
3D, and considerable theoretical uncertainty of the under- 
lying fluid physics. With the above caveat in mind, we will 
use as a flrst approximation the toy surface density evolution 
model introduced in ii3.2l to estimate the thermal emission 
from an accreting PTA source. 

In Figure [S] we show the thermal spectrum for a Mg = 
1, M2/M1 = 1/4 PTA source, calculated from the circumbi- 
nary gas surface density profiles in Figure|4l We have plotted 
(i) the circumbinary disc with rh — 3 inside the radius where 
Q = 1 (dotted, left hump); (ii) an accretion disc around the 
secondary SMBH (dotted, right hump) fueled by leakage into 
the cavity and truncated at the Hill radius, for which we use 
Rh ~ 0.5r7^'''a; and (iii) the combined emission of the two 
discs (solid thick line). For comparison, we also show the 
spectrum for an Eddington-limited thin disc (dashed lines) 
around a solitary SMBH with the same mass as the binary. 
For simplicity, we have assumed that all of the gas leaked 
into the cavity fuels a circum-secondary disc. We have plot- 
ted spectra when the source has a binary orbital period of 
P = 1 yr and when P = 0.1 yr. 

The result is what would be expected intuitively. The 



infrared and optical flux, which is produced almost exclu- 
sively in the circumbinary disc, does not vary greatly from 
what is expected from a standard thin disc. However, the 
flux drops precipitously below wavelengths of A ;$ SOOOA 
{u > 10^^ Hz in the figure). This is in stark contrast to most 
unobscured quasars thought to be powered by ~ W^~^ Mq 
SMBHs, which have their brightest emission in the rest- 
frame near-UV near their Lyman-a line. The bolometric lu- 
minosity of the accreting PTA source is roughly ~ O.OSLEdd 
(i.e., L/LEdd ~ 10~^m) for P = 1 yr, and ~ lO^'^Z/Edd for 
P = 0.1 yr. The optical and infrared emission is dominated 
by the circumbinary disc, whereas the UV and X-rays are 
produced by circum-secondary accretion fueled by leakage of 
circumbinary gas into the cavity. As the binary evolves to- 
ward shorter periods, the circum-secondary disc is depleted 
— the viscous time at the Hill radius is typically a few hun- 
dred years, shorter than the time to binary merger — and 
as a result, less gas is able to leak into the cavity, decreasing 
the UV and X-ray emission. The degree to which the UV 
and X-ray emission is suppressed depends on the model pa- 
rameters (in particular /leak) and on the binary period. Note 
that the system may still be luminous in hard X-rays due 
to inverse Compton scattering by a coronal electron plasma 
(|Sesanaet al.ll201ll ). 

We also note that the downturn in the near-UV flux at 
A ^ 300 nm could help distinguish PTA sources from single- 
SMBH AGN. This feature will be observable in the optical 
if the source redshift is high; e.g., at z = 1 it will be in the 
Vband. Hence, even in the optical, this source will have an 
unusual color: it will appear fainter in the U and B bands 
than a typical AGN. The downturn could be distinguished 
from reddening due to dust obscuration through the devia- 
tion from the power-law spectral shape of dust reddening. 

We propose that once an individually resolved PTA 
source is detected and its error box determined, search- 
ing for AGN with weak UV emission lines (e.g., Ly a) 
and/or weak soft X-ray emission is a promising method 
to narrow the field of interlopers. AGN whose soft X-ray 
fluxes are weaker by more than a factor of 10 compared to 
the average have indeed been detected, and are estimated 
to constit ute at most ~ 1% o f the ge neral AGN popula- 
tion (e.g.. iBrandt. Laor fc W ills 2000 !: iLeighlv et al.l [2OO7I : 
iGibson. Brandt fc Schneideill2008. : iWu et al.ll201ll . and refs. 
therein). There have als o been observations of quasa rs with 
exceptionally weak lines (|Diamond-Stanic et al.ll2009l ): these 
objects have infrared and optical emission consistent with 
those of typical luminous A GN, and also tend to be X-ray 
weak (|Shemmer et al.ll2009l ). That X-ray weak AGN are so 
rare suggests that it will be possible to narrow the num- 
ber of interlopers in a typical PTA error box by a factor of 
~ 100, i.e. either to a handful of objects, or yielding a unique 
EM counterpart candidate. It is possible, furthermore, that 
some of these rare X-ray weak AGN are in fact the SMBH 
binaries that PTAs will be detecting. 

Our results also strongly suggest that AGN counter- 
parts to PTA sources should draw from optically selected 
surveys, as their nature mak es them likely to b e missed by 
X-ray searches (see, however, ISesana et al.ll201ll , who inves- 
tigate the possible X-ray searches of PTA source binaries 
that have not yet decoupled). The Large Synoptic Survey 
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Telescopiij should be able to detect all of the optically lu- 
minous AGN in the PTA error box within z ~ 1. It may be 
possible to follow up candidates individually, but comparing 
the optical data with that of wide-field X-ray surveys such as 
MAXE, eROSITy^ or Wide Field X-ray Telescop^ would 
greatly facilitate the multi-wavelength search for counter- 
part candidates inside the error box. 

Additional follow-up studies of candidates may further 
corroborate the identification of a counterpart. For exam- 
ple, the gas that leaks radially into the cavity can shock- 
heat the outer edge of the circum-secondary (or circum- 
primary) disc and produce hot spots. The viscously dis- 
sipated luminosity of a circum-secondary disc is roughly 
LdiBc2 ^ {l/2)GM2M2/Risco,2, where Ma < /ioakM(i?A) 
is the mass supply rate of the circum-secondary disc and 
-Risco,2 is the radius of innermost stable circular orbit 
around the secondary. The time-averaged power per unit 
mass of the hot spots is limited by the amount of kinetic 
energy the flow can deposit at the outer edge of the circum- 
secondary disc, i.e. Lhot ^ GlVhA-h/RH- It follows directly 
that the time-averaged ratio between between the hot-spots 
and the intrinsic luminosity of the circum-secondary disc is 



equation 15) 



Lhot ^ flisco,2 



-l/3^ISCO,2 ^ / 



La 



Rh 



\GM/c^ 



(43) 



In other words, the time-averaged power of a hot spot is of 
order ^ 1 — 10% of the circum-secondary disc luminosity for 
resolved PTA sources. In principle, the luminosity of any sin- 
gle flare could be much greater. Because streaming into the 
cavity is expected to be mod ulated quasi-periodically by the 



binary's orbital period (e.g., Havasaki. Mineshige fc Sudoul 



I2OO7I : iMacFadven fc Milosavlievid l2008l l. EM counterparts 



of resolved PTA sources may be characterized by periodic 
UV flares. 

In the same vein, if the orbital plane lies close to 
the line of sight, the UV lines would display strong 
periodic Doppler shifting with respect to the optical 
emission, modulated at the binary's orbital period (e.g., 
iHalpern fc Filippenkd Il98g) . Thus, monitoring candidate 
counterparts for periodic or quasi-periodic variability on or- 
bital timescales may prove a fru itful route for identification 
iHaiman. Kocsis fc Men ou \2QQ9) . As a proof of this concept, 
we note that iBoroson fc Laucr ((2009|) recently reported a 
candidate SMBH binary, with two sets of broad emission 
lines separated by 3,500km s~^, with inferred component 
masses of Mi = lO*'® Mq and M2 = lO'''^ Mq. The binary 
interpretation, however, could be ruled out by the lack of 
any change in t he velocity offset betw een two spectra taken 
« 1 year apart (jChornock et al.ll2009l ). 

Lastly, we consider the scenario of IChang et al.l (|201(J ). 
in which the circum-primary disc brightens prior to merger 
due to tidal excitation by the shrinking binary. The power 
generated by this process can be approximated as (see their 



'' http://www.lsst.org/lsst 

* http://maxi.riken.jp/top/ 

^ http://www.mpc.mpg.de/crosita/ 

^'^ http://www.wfxt.eu/homc/Overview.html 
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where Min is t he mass of the circu m-primary disc. Extrap- 
olation of the IChang et al.l (|20ld ) results to binaries with 
mass M ~ 10^ Mq (their calculations only considered bina- 
ries up to M = lO^Mg) suggests a value of Min ~ IOOMq. 
Similar values are obtained by estimating the disc mass 
that can be fueled by gas leaking into the cavity with 
/leak ~ 0.1, when the time to merger is comparable to the 
viscous timescale at the outer edge of the circum-primary 
disc. 

Equation 1441 suggests that the power produced by tidal 
excitation of the circum-primary disc is negligible compared 
to the thermal disc emission if the binary period is P ~ 1 yr. 
However, for sources with P ~ 0.1 yr, the tidally excited 
emission would rival the bolometric output of the thermal 
emission. The tidal component, which would have a peak 
frequency in the UV and soft X-rays, will brighten dra- 
matically prior to merger on timescales of several years to 
decades. Even though P ~ 0.1 yr sources are predicted to 
be much rarer and also much more difficult to resolve indi- 
vidually with PTAs, they present tantalizing possibilities for 
observing EM signatures that are directly related to binary 
coalescence. 



4 CONCLUSIONS 

In this paper we considered the possibility that individually 
resolved PTA sources — SMBH binaries with M ~ lO'^M©, 
M2/M1 ~ 1/4, P ~ 0.1 - 1 yr and z <, 1.5 — may be 
identified with EM observations if they reside in gas-rich 
environments. Multi-wavelength observations of such sys- 
tems would allow for studies of an AGN in a system that 
is known to harbour a compact SMBH binary, thus pro- 
viding a unique window into gas accretion in a rapidly 
time-varying grav i tation al potential. If, as suggested by 
ICorbin fc CornishI ((20101), PTAs can constrain the luminos- 
ity distance to individually resolved sources, these SMBHs 
can be used as "standard sirens" to measure the cosmic ex- 
pansion history. Interestingly, the predicted redshift distri- 
bution of these sources lies between 0.1 ^ z ^ 1.5, a range 
comparable to that of the deepest Type-la supernovae sur- 
veys and where LISA detection rat es are expected to be low 
ISesana. Volonteri fc HaardjIJOOTJ ). 

Our findings can be summarized as follows: 

• The number of interloping massive halos and AGN that 
may be confused with the PTA source is typically Ng ~ 10* 
for a 10^ M0 binary if the contributions to the signal from 
individual pulsars are not identifiable in the PTA data. 

• In the more optimistic case, the pulsar term can be 
constrained and utilized to better determine the sky location 
of the source as well as constrain its redshift and mass. The 
number of interlopers can then drop to Ng ~ 10 — 100 at 
z ~ 0.5, and perhaps to N < 1 for close sources at redshifts 
as close as 2: ^ 0.2. 

• By considering the orbital evolution history of an ac- 
creting PTA source, first by tidal interactions with the cir- 
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Figure 5. We plot the spectral energy distribution from the thin, 
viscous circumbinary accretion disc from Figure |4] (red curve; 
M = lO^AfQ, P = 1 yr), when the source has a period (a) 
P = 1 yr and (b) P = 0.1 yr. We assume that some of the 
gas at the cavity/disc boundary is able to accrete into the cavity 
at a suppressed rate and fuel a small accretion disc around the 
secondary SMBH (see text for details). The spectra from individ- 
ual discs are shown in dotted lines, with the circumbinary cavity 
emitting at lower frequencies. The combined spectrum from the 
circumbinary and circum-secondary accretion discs is shown in 
solid thick lines. We plot for reference the model AGN spectrum 
(thin dashed) for an Eddington-limited thin disc around a single 
SMBH with the same mass as the binary. The spectra for the PTA 
source are UV and X-ray weak compared to a standard thin disc 
model around a single SMBH. However, the infrared and optical 
emission, mostly produced in the circumbinary disc, is similar to 
what would be expected for a single-SMBH disc. 



cumbinary gas and then by GW emission, we showed that 
a gaseous accretion disc around the source can be expected 
to be gas-poor both inside and immediately outside its or- 
bit. They would thus have optical and infrared luminosities 
comparable with typical quasars, while exhibiting low ther- 
mal X-ray luminosities and weak UV emission lines. The 
downturn in flux below 300 nm could be discerned by opti- 
cal observations if the source redshift is z > 1. The leakage of 
circumbinary gas into the cavity may shock the circumpri- 
mary (-secondary) disc and cause substantial quasiperiodic 
fluctuations in the UV and X-ray. Searching for AGN in 
the PTA error box with one or more of these atypical char- 
acteristics could lead to the identification of a single EM 
counterpart. Further monitoring candidate counterparts for 
periodicity and other theoretically predicted pre-coalescence 
signatures may also aid identification. 

• The above emission features are general for AGN pow- 
ered by thin accretion discs around compact SMBH binaries. 
They could thus be used to discover such systems even in the 
absence of a GW detection. This is a particularly interesting 
possibility, as a small fraction (^ 1%) of optically luminous 
AGN is known to have low X-ray luminosities and unusually 
weak emission lines. Upcoming wide-field surveys in the op- 



tical and X-rays should together discover many more such 
AGN, as well as observe them over temporal separations 
of weeks to months. Such cadential, multi-wavelength data 
may lead to the first observations of compact SMBH binaries 
with sub-parsec separations. 

We have assumed that a PTA source would appear as 
luminous AGN, with the gaseous fuel perhaps being sup- 
plied by the preceding merger of the binary's host galaxies. 
The degree to which galaxy mergers dictate AGN activity 
remains an open question, and the link between SMBH bi- 
narity and AGN activity is even less certain. It is possible 
that many SMBH binaries that are individually detected by 
PTAs will have no EM counterpart at all. 

Our results were calculated using a simple semi-analytic 
accretion disc model, a central assumption of which is that 
the binary's tidal torques are able to open a central cavity 
in the disc. In the radiation-dominated regions of interest, 
strong horizontal advective fluxes or vertical thickening of 
the disc may act to close such a cavity and wipe out the 
features we predict. The features would also not be present 
if the disc and binary's orbits do not lie on the same plane, 
as i n the binary model of t he variable BL Lac object OJ 
287 (jLehto fc Valtonenlll996l ). Absorption and reprocessing 
by the binary's host galaxy may also act to mask or mimic 
the intrinsic thermal AGN emission we have modeled. 
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APPENDIX: GREEN'S-FUNCTION SOLUTION FOR THE THIN-DISC EQUATION WITH MOVING 
INNER BOUNDARIES 

Equation [35] states that given an initial surface density profile Einit(-R')i ^he general solution E(_R, i) to the Keplerian thin- 
disc equation (equation I34|) can be written in the form E(i?) = J G{R,R') Einit(fl') dR' . The Green's function G may be 
thought of as a transform kernel that solves the thin-disc equation for any choice of the time i, inner boundary Rx > 0, the 
viscosity power-law n = lni//ln_R < 2 and the "leakage parameter" /icak (defined in equation I36p that quantifies the mass 
flow rate across R\. Below, we consider a case where where the boundary R\ moves inward as an explicit function of time. 
For simplicity, we assume n and /loak to be constants. 

Our astrophysical problem has two convenient properties that aid our mathematical analysis: first, that the gas leaking 
into the cavity does so in nearly radial orbits means it is effectively decoupled from the circumbinary disc and the matter 
inside R < Rx can be ignored; second, the disc is empty inside the initial value of Rx and so it is only necessary to integrate 
over R' > Rxfi = Rx{t = 0). Our present approach will not work, for example, for a general boundary that moves outward 
with time. 

Suppose that Q{R,R';t*) is a weighted sum of the Green's function G that consists of solutions evaluated at different 
values of Rx{t) over the duration < t < t* . By the principle of superposition, Q must also be a solution of the thin-disc 
solution that satisfies 



T.{R,t') 



g{R,R';t') Sinit(i?') dR'. 



(45) 



At any t and Rx, G must satisfy dQ jdt = dG/dt and dQ/dRx = dG/dRx- Further, if Rx is an explicit function of time, it 
must be true that the full time dependence of Q is described by dQ/dt — dG/dt -\- {dG / dRx){dRx / dt) . It then follows that 



E(i?,i + At) 



GiR, R'; t, Rx) + ^^l^i^jliMAt + O(Af^) + 



Einit(i?') dR' 



GiR,R';t,Rx)+ ^^^^'^':^'^'^ At+ ^^^\^y'^''^ ^At+ . 



dt 



dRx 



dt 



Einit(-R ) dR . 



From equations 1351 and 1461 we may write the derivative dT,/dt through its definition: 



E(i?,t) 



dG{R, R';t, Rx) dG{R, R'; t, Rx) dRx 



dt 



dRx 



dt 



Einit(i?') dR!. 



Direct integration gives the expression 

E(ii:,f*) = Ei„it(E') + 

JO JRxfi 

Our "master" Green's function G is therefore 

G{R,R!;t') ^S{R~R') + 



dG{R, R';t, Rx)_ dG(R, R'; t, Rx) dRx 



dt 



dRx 



dt 



Einit(i?') dR' dt. 



dG{R, R'; t, Rx)_ dG{R, R'; t, Rx) dRx 



dt 



dRx 



dt 



dt. 



(46) 



(47) 



(48) 



(49) 



The Dirac (5-function can be awkward to implement in a numerical integration scheme. We rewrite it in terms of the known 
Green's function with a fixed boundary by using the fact that any Green's function evaluated at i = is the (5-function, i.e.. 



GiR,R';t'',R*x) = S{R-R 



')+ f 
Jo 



di 



dt, 



where R\ — Rx{t*). 

We write our solution as 



G{R,R';t')^G{R,R',t'-Rl) + 



t* r 



oCryR, R ; i, Rx 



dt 



dG{R, R';t, R*x)_ dG(R, R'; t, Rx) dRx 



dt 



dRx 



dt 



dt. 



(50) 



(51) 



We see that if the boundary is stationary, the second term above vanishes an d G = G. Note the similarity of the math- 
ematical form of our solution to DuHamel's theorem fe.g.. lCarslaw fc Jaegeilll959l ) for time-dependent boundary conditions. 
This is not surprising, given that both are based on the superposition principle; equation [51] is simply a weighted sum of 
instantaneous Green's functions and not, strictly speaking, a new type of solution. For our thin accretion disc problem, the 
function G{R,R',t) is explicitly known and easily tabulated, given a specific combination of the parameters n and /icak and 
the boundary evolution Rx{t). 

In the context of this paper, the initial condition is the point when the binary's orbital decay becomes GW-driven, and t* 
is the time when the system is observed as an individually resolved PTA source and Rx{t) > Rx describes the (GW-driven) 
evolutionary history of the gap-opening radius prior to t* . 
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